!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!>\brief
!! Error computations for the Stokes/Navier-Stokes problem.
!!
!! \n
!!
!! All the subroutines of this module assume that an analytic solution
!! is available. Errors are computed with numerical quadrature, and
!! it's advisable to use more accurate quadrature formulas than those
!! used in finite element method. The extra accuracy is controlled by
!! the input argument \c extra_deg used by all the error computation
!! subroutines.
!<----------------------------------------------------------------------
module mod_error_norms

!-----------------------------------------------------------------------

 !$ use omp_lib

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_grid, only: &
   mod_grid_initialized, &
   t_grid, affmap

 use mod_base, only: &
   mod_base_initialized, &
   t_base, &
   new_base => base

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_error_norms_constructor, &
   mod_error_norms_destructor,  &
   mod_error_norms_initialized, &
   velocity_err, &
   pressure_err

 private

!-----------------------------------------------------------------------

! Module types and parameters

! Module variables

 ! public members
 logical, protected ::               &
   mod_error_norms_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_error_norms'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_error_norms_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
       (mod_kinds_initialized.eqv..false.)    .or. &
       (mod_grid_initialized.eqv..false.)     .or. &
       (mod_base_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_error_norms_initialized.eqv..true.) then
   endif
   !----------------------------------------------

   mod_error_norms_initialized = .true.
 end subroutine mod_error_norms_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_error_norms_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_error_norms_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_error_norms_initialized = .false.
 end subroutine mod_error_norms_destructor

!-----------------------------------------------------------------------

 !> Error computation for the velocity variable. The following errors
 !! are computed:
 !! <ul>
 !!  <li>\c e_l2
 !! \f$=\left( \sum_K 
 !!   \|\underline{u} - \underline{u}_h\|^2 \right)^{\frac{1}{2}} \f$
 !!  <li>\c e_h1
 !! \f$=\left( \sum_K
 !!   \|\nabla\underline{u} - \nabla\underline{u}_h\|^2 + e_{L^2}^2
 !!           \right)^{\frac{1}{2}} \f$
 !!  <li>\c e_div
 !! \f$=\left( \sum_K \|\nabla\cdot\underline{u}_h \|^2
 !! \right)^{\frac{1}{2}} \f$
 !! </ul>
 !! where \f$K\f$ are the mesh elements.
 !!
 !! \note For consistency with the other subroutines in this module
 !! the format of uuuh is the same as for a discontinuous method, with
 !! element based degrees of freedom corresponding to the first index
 !! and spatial dimension corresponding to the second index (i.e. as a
 !! collection of scalar discontinuous fields).
 !!
 !! \bug The interface of \c uuu should be defined using the abstract
 !! interfaces given in \c mod_testcases.
 !<
 subroutine velocity_err(base_in,grid,uuuh,uuu,guuu,extra_deg, &
   e_l2,e_h1,e_div)
  type(t_base), intent(in) :: base_in
  type(t_grid), target, intent(in) :: grid
  real(wp), intent(in) :: uuuh(:,:)
  interface
    pure function uuu(x) result(u)
     import :: wp
     real(wp), intent(in) :: x(:,:)
     real(wp) :: u(size(x,1),size(x,2))
    end function uuu
  end interface
  interface
    pure function guuu(x) result(gu)
     import :: wp
     real(wp), intent(in) :: x(:,:)
     real(wp) :: gu(size(x,1),size(x,1),size(x,2))
    end function guuu
  end interface
  integer, intent(in) :: extra_deg
  real(wp), intent(out) :: e_l2, e_h1, e_div

  integer :: ie, i, j
  real(wp) :: l2, h1
  real(wp), allocatable :: xg(:,:), wg(:), uuue(:,:), grad_uuue(:,:,:), &
    uuuhk(:), uuuhe(:), guuuhe(:,:), div(:)
  type(t_base) :: base ! basis with higher accuracy
  character(len=100) :: message(3)
  character(len=*), parameter :: &
    this_sub_name = 'velocity_err'

   ! increase the order of the quadrature formula
   base = new_base( base_in%me , p_s=base_in%p_s , &
                    deg=base_in%deg+extra_deg )

   allocate( xg(base%me%d,base%m), wg(base%m),                      &
     uuue(base%me%d,base%m), grad_uuue(base%me%d,base%me%d,base%m), &
     uuuhk(base%pk), uuuhe(base%m), guuuhe(base%me%d,base%m),       &
     div(base%m) )

   e_l2  = 0.0_wp
   e_h1  = 0.0_wp
   e_div = 0.0_wp
   elem_loop: do ie=1,grid%ne
     xg = affmap(grid%e(ie),base%xig)
     wg = grid%e(ie)%det_b * base%wg

     ! analytic solution
     uuue = uuu(xg)
     grad_uuue = guuu(xg)

     ! numerical solution and error computation
     l2 = 0.0_wp
     h1 = 0.0_wp
     div = 0.0_wp
     do i=1,base%me%d
       ! get the numerical solution
       uuuhk = uuuh((ie-1)*base%pk+1 : ie*base%pk,i)

       ! L2 error
       uuuhe = matmul(uuuhk,base%p)
       l2 = l2 + sum(wg*(uuue(i,:)-uuuhe)**2)
       
       ! H1 error
       do j=1,base%me%d ! compute \nabla u_i
         guuuhe(j,:) =  matmul(uuuhk,base%gradp(j,:,:)) ! du_i/dx_j
       enddo
       guuuhe = matmul(transpose(grid%e(ie)%bi),guuuhe)
       do j=1,base%me%d
         h1 = h1 + sum(wg*(grad_uuue(i,j,:)-guuuhe(j,:))**2)
       enddo

       ! div error
       div = div + guuuhe(i,:)

     enddo
     e_l2  = e_l2  + l2
     e_h1  = e_h1  + l2 + h1
     e_div = e_div + sum(wg*(div**2))

   enddo elem_loop

   deallocate( xg, wg, uuue, grad_uuue,  &
               uuuhk, uuuhe, guuuhe, div )

   if( (e_l2.lt.0.0_wp).or.(e_h1.lt.0.0_wp).or.(e_div.lt.0.0_wp) ) then
     write(message(1),'(a,e10.3,a,e10.3,a,e10.3,a)') &
       'One error is negative: e_l2**2 = ',e_l2,', e_h1**2 = ',e_h1, &
       ', e_div**2 = ',e_div,';'
     write(message(2),'(a)') &
       '  A possible reason is the use of a quadrature formula with'
     write(message(3),'(a)') &
       '  negative weights. Error set to zero.'
     e_l2  = max(e_l2,0.0_wp)
     e_h1  = max(e_h1,0.0_wp)
     e_div = max(e_div,0.0_wp)
     call warning(this_sub_name,this_mod_name,message)
   endif

   e_l2  = sqrt(e_l2)
   e_h1  = sqrt(e_h1)
   e_div = sqrt(e_div)

 end subroutine velocity_err

!-----------------------------------------------------------------------

 !> Error computation for the pressure variable. The following errors
 !! are computed:
 !! <ul>
 !!  <li>\c e_l2
 !! \f$=\left( \sum_K \|p - p_h\|^2 \right)^{\frac{1}{2}} \f$
 !!  <li>\c e_h1
 !! \f$=\left( \sum_K
 !!   \|\nabla p - \nabla p_h\|^2 + e_{L^2}^2
 !!           \right)^{\frac{1}{2}} \f$
 !! </ul>
 !! where \f$K\f$ are the mesh elements.
 !!
 !! \note For consistency with the other subroutines in this module
 !! the format of ppph is the same as for a discontinuous method.
 !!
 !! \bug The interface of \c ppp should be defined using the abstract
 !! interfaces given in \c mod_testcases.
 !<
 subroutine pressure_err(base_in,grid,ppph,ppp,gppp,extra_deg,e_l2,e_h1)
  type(t_grid), target, intent(in) :: grid
  type(t_base), intent(in) :: base_in
  real(wp), intent(in) :: ppph(:)
  interface
    pure function ppp(x) result(p)
     import :: wp
     real(wp), intent(in) :: x(:,:)
     real(wp) :: p(size(x,2))
    end function ppp
  end interface
  interface
    pure function gppp(x) result(gp)
     import :: wp
     real(wp), intent(in) :: x(:,:)
     real(wp) :: gp(size(x,1),size(x,2))
    end function gppp
  end interface
  integer, intent(in) :: extra_deg
  real(wp), intent(out) :: e_l2, e_h1

  integer :: ie, i
  real(wp) :: l2, h1
  real(wp), allocatable :: xg(:,:), wg(:), pppe(:), &
    ppphk(:), ppphe(:), gpppe(:,:), gppphe(:,:)
  type(t_base) :: base ! basis with higher accuracy
  character(len=100) :: message(3)
  character(len=*), parameter :: &
    this_sub_name = 'pressure_err'

   ! increase the order of the quadrature formula
   base = new_base( base_in%me , p_s=base_in%p_s , &
                    deg=base_in%deg+extra_deg )

   allocate( xg(base%me%d,base%m), wg(base%m),         &
     pppe(base%m), ppphk(base%pk), ppphe(base%m),      &
     gpppe(base%me%d,base%m), gppphe(base%me%d,base%m) )

   e_l2 = 0.0_wp
   e_h1 = 0.0_wp
   elem_loop: do ie=1,grid%ne
     xg = affmap(grid%e(ie),base%xig)
     wg = grid%e(ie)%det_b * base%wg

     ! analytic solution
     pppe  = ppp(xg)
     gpppe = gppp(xg)

     ! numerical solution
     ppphk = ppph((ie-1)*base%pk+1 : ie*base%pk)
     ppphe = matmul(ppphk,base%p)

     ! L2 error
     l2 = sum(wg*(pppe-ppphe)**2)

     ! H1 error
     do i=1,base%me%d ! compute \nabla p
       gppphe(i,:) =  matmul(ppphk,base%gradp(i,:,:))
     enddo
     gppphe = matmul(transpose(grid%e(ie)%bi),gppphe)
     !      int           point-wise error
     h1 = sum( wg * sum((gpppe-gppphe)**2,1) )

     e_l2 = e_l2 + l2
     e_h1 = e_h1 + l2 + h1

   enddo elem_loop

   deallocate( xg, wg, pppe, ppphk, ppphe, gpppe, gppphe )

   if( (e_l2.lt.0.0_wp).or.(e_h1.lt.0.0_wp) ) then
     write(message(1),'(a,e10.3,a,e10.3,a)') &
       'One error is negative: e_l2**2 = ',e_l2,', e_h1**2 = ',e_h1,';'
     write(message(2),'(a)') &
       '  A possible reason is the use of a quadrature formula with'
     write(message(3),'(a)') &
       '  negative weights. Error set to zero.'
     e_l2 = max(e_l2,0.0_wp)
     e_h1 = max(e_h1,0.0_wp)
     call warning(this_sub_name,this_mod_name,message)
   endif

   e_l2  = sqrt(e_l2)
   e_h1  = sqrt(e_h1)

 end subroutine pressure_err

!-----------------------------------------------------------------------

end module mod_error_norms

